Viscosity and mutual diffusion in strongly asymmetric binary ionic mixtures 
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We present molecular dynamics simulation results for the viscosity and mutual diffusion constant 
of a strongly asymmetric binary ionic mixture (B1M). We compare the results with available theo- 
retical models previously tested for much smaller asymmetries. For the case of viscosity we propose 
a new predictive framework based on the linear mixing rule, while for mutual diffusion we discuss 
some consistency problems of widely used Boltzmann equation based models. 
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The important advancements that occurred in the last decade in the experimental techniques involving high-power 
lasers have generated a renewed interest in the properties of dense plasmas in general and their transport properties 
in particular 0, 0, 0] • The experimental capabilities currently available and the ones that are expected to become 
available in the near future Q promise to further advance the field of inertial confinement fusion (ICF) as well as shed 
new light on long-standing astrophysics problems. Many times such experiments either probe, or their outcome is 
O |' strongly dependent on, the behavior of plasma mixtures with various degrees of charge and mass asymmetries of the 
components. Such mixtures occur for example in ICF due to the instability (e.g. Richtmyer-Meshkov or Rayleigh- 
Taylor) driven mixing of the heavy elements that make up the enclosing shell and the much lighter fuel. In this case 
the stability of the initial interfaces, nature of the ensuing flows and degree of fuel contamination are crucially linked 
to such mixture properties as shear viscosity and mutuaf diffusion 0,10. In the present paper we calculate these 
properties using molecular dynamics simulations for a simple but relevant plasma model: the binary ionic mixture 
(BIM) which is a generalization of the one-component plasma (OCP) [13 . We study a rather extreme, ICF- 

inspired system 0, D + — Au 39+ , which displays roughly two orders of magnitude charge and mass asymmetry, but 
the results should apply to other plasma mixtures as well, e.g. of astrophysics significance [Tlj . where large charge 
>^Y and mass asymmetries are present. We compare the MD results with available theoretical models previously tested 
for much smaller asymmetries. For the case of viscosity we propose a new predictive framework based on the linear 
mixing rule, while for mutual diffusion we discuss some consistency problems of Boltzmann equation based models. 

The BIM model consists of a mixture of Ni point ions of charge q± = Z x e and mass Mi and N 2 point ions of charge 
12 = Z 2 e and mass M 2 embedded in a uniform, rigid, neutralizing electronic background. We denote the number 
concentrations by x a — N a /N, N = Ni + N 2 , a = 1,2 and number densities by p a — N a /V, p = pi + p 2 , where 
V — L 3 is the volume. (Z) — X\Z\ + x 2 Z 2 is the average charge and pi = p\Z\ + p 2 Z 2 the electronic number density. 
As usual the mean inter-electronic and inter- ionic distances are defined by ai = (3/ 'An pi) 3 and a — (i/Airp)^ = ai(Z)^ , 
while the electronic and ionic coupling parameters are: 
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The thermodynamics of the system is fully determined by one coupling constant, e.g. T, and concentration x±. As it 
is the case for the OCP, the validity range of the BIM model is such that the Fermi temperature of the electrons is 
Tp T and the densities are high enough so r s <C 1, r s — ai/ao, ciq - Bohr radius, corresponding to a completely 
" ' degenerate and rigid electronic background. 

The thermodynamics of the BIM has been thoroughly studied and is known to be very well described by the linear 
mixing rule [§1 u3. IT3L IT^. For moderate charge asymmetries an OCP-based "one-fluid theory" is also a reasonable 
approximation with an effective charge = suggested by the ion-sphere model [l3 |. The relative 

success of this "one-fluid" representation has lead Clerouin et al. to propose that the shear viscosity can also be 
predicted in terms of the equivalent OCP, as already tested for thermal transport ^(|. (Although not explicitly 
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stated in , further assumptions need to be made for such a prediction - see below) . This "one-fluid" approach was 
shown to be suitable for calculating the BIM viscosity at charge and mass asymmetries of order « 10. Before testing 
this idea on much larger asymmetries, « 100, we note that, surprisingly, the viscosity of the OCP itself does not 
appear to be very accurately known. For intermediate and strong couplings, 1 < T < 100, Bernu and Vieillefosse 
fl7] have proposed an interpolation formula based on three MD simulation results obtained with systems of 128 - 
250 particles, while at even stronger couplings Ichimaru and Tanaka have introduced a generalized viscoelastic theory 
El. In the authors propose a different relation based on the kinetic theory of Wallenborn and Baus E3i that 
extends to the weak coupling regime. However, the disagreement between these two approaches (Refs. E3 and @|) is 
significant in the regime that they both cover, T > 1, particularly at intermediate T's. Unfortunately it is difficult to 
ascertain the reliability of these predictions given both the limited simulation results available and the small system 
size used, which limits the accuracy of the results. 

To settle this question we performed extensive microcanonical MD simulations of the OCP with much larger system 
sizes - 1372 particles, and a wide range of coupling constants, 0.05 < T < 100. The Coulomb interactions were handled 
using the Ewald summation technique with conducting boundary conditions. The calculation of the shear viscosity r\ 
was done using the Green-Kubo relation: 
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As shown by Bernu and Vieillefosse 17], and more recently in the context of Yukawa plasmas by Salin and Caillol 
|20|. the evaluation of the pressure tensor a requires an Ewald- type summation for its interaction part, a 1 : 
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where a and b denote the Cartesian coordinates. The Ewald parameter a [2lJ was chosen such that the real space 
sums, e.g. a^ r \ can be calculated with the usual minimum-image convention, as shown above. The duration of the 
runs was 10 3 — lO 4 ^" 1 , where lo 2 = Airpe 2 /M is the plasma frequency [Z = 1). The natural unit for the viscosity of 
the OCP is T] Q = pMa 2 uj p . 

Our simulation results are shown in Fig. ^ together with the interpolation formula of Bernu and Vieillefosse (for 
r > 1) and the relation proposed in based on OCP kinetic theory. The errors, estimated using a standard block 
analysis (2^, are between about 8% at intermediate and large T and 25% at the lowest V. The present viscosity results 
largely agree with those of [17|], but suggest that in the intermediate coupling range the OCP viscosity is significantly 
higher than previously predicted, in agreement with |23|. The well known viscosity minimum appears to be around 
r ~ 21, with rj/r/o ~ 0.084. We are not aware of other simulations for weakly coupled plasmas, T < 1, but our viscosity 
results in this regime are in qualitative agreement with the kinetic theory of Wallenborn and Baus, although somewhat 

lower. Since for these conditions the screened, Debye-Hiickel potential e 2 exp(—r/\o)/r, \d = [fcsT/47re 2 p] 2 , should 
be an appropriate representation of the effective inter-ionic interaction and a/Xjj oc T^, i.e. the interaction can 
be fairly short-ranged even for T significantly smaller than 1, a good viscosity estimate should be provided by the 

Chapman- Enskog theory, 77 = 5kgT/8fl^ [24[ . As shown in Fig. ^(see further below for how the collision integrals 

(2) 

ill, are calculated) this is indeed the case for 0.05 < T < 1. We choose to fit all the data points with the relation: 



= AT- 2 + Br~ s + CT 



(11) 



which captures rather well the behavior of the OCP viscosity in the wide range of couplings simulated. The best 
parameters are: A = 0.482, B = 0.629, C = 1.88 x 10~ 3 , and s = 0.878. 
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We know tackle the question of the BIM viscosity, in particular when the charge and mass asymmetries are very 
large. It should be noted that the system that we focus on, D + — Au 39+ , is not a simple BIM per se, as gold {Au) 
is only partially ionized and the effect of the remaining electrons may be important under certain thermodynamic 
conditions. However, the BIM approach is still relevant provided the densities and temperatures are such that the 
distance of closest approach between ions is larger that the radius of the remaining ion cores. 

To elucidate the effect on viscosity of mixing plasmas with such large differences in charge and mass as D + and 
Au 39+ it is convenient to adopt the procedure of Ref. 0, where the coupling constant T is kept fixed and the 
concentration of the two species is varied. We set T = 0.05, which corresponds to I\ = 0.05 for pure deuterium 
(xau = 0), i.e. a weakly coupled plasma, and to I\ ~ 76 for pure gold (xa u = 1), i-e. a strongly coupled plasma. As 
before, we perform microcanonical simulations with a system of 1372 particles, at a number of different concentrations 
xau- Due to the strong charge and mass asymmetries, exceedingly long run times are necessary for both equilibration 
and data accumulation to calculate the viscosity using Eq. [Hwith an accuracy of 20 — 25%. A good measure of the 
large "size" difference between the ions is provided for example by the pair correlation functions, which we show in 
Fig. Elfor x A u = 0.5. 

As noted in |2j the viscosity drops steeply upon mixing highly charged, heavy ions in a weakly coupled plasma, 
see Fig. [3J This effect can be understood qualitatively in the framework of a "one-fluid" theory, which we describe 
below. The coupling constant of the equivalent OCP is T e ff = ^^eff' wnere %eff = follows from the 

ion-sphere model |l5j . In this approximation the thermodynamics is fully determined by ^ e ff, but the calculation of 
the viscosity requires some additional arguments. For example, a reasonable unit for the viscosity of this system may 
be assumed to be 77™ = p(M)a 2 Lu pm , (M) — X\M% +X2M2, where ujp m = ui'j(Z)' 2 / (M) is the "hydrodynamic" plasma 
mixture frequency. (The use of the so-called "kinetic" mixture frequency ]25| leaves the results largely unchanged) . 
The mixture viscosity in these units is then postulated to be given by the scaled effective-OCP viscosity: 

v(T,x) = Vocp(T e ff) , 12 , 

We show in Fig. [3] the results of such calculations for the system that we study, using both the OCP viscosity of 
Ref. 0], and our prediction Eq. ^2 The qualitative dependence of the MD results on xaui arising from the opposite 
behaviors of rjocpi^eff) and 77™ (x), is reproduced correctly, but the quantitative disagreement is also very significant, 
particularly for small and moderate Au concentrations. The use of the more accurate OCP viscosity Eq. ^]does not 
fully alleviate this problem. We conclude that, not surprisingly, the accuracy of the "one- fluid" model is diminished 
for extreme asymmetries. 

It is clear that this limitation can only be overcome by taking into account, either explicitly or implicitly, the 
mixture asymmetry. A direct calculation along the lines of the kinetic theory of Wallenborn and Baus [jjj] has 
been used for example to determine the BIM mutual diffusion constant at small asymmetries |2(|. However, the 
calculation of the viscosity is even more complex and given the limitations of the theory even for the OCP its success 
for BIM quantitative predictions is rather doubtful. We turn therefore to a more indirect method, which we outline 
below. First, we recall an interesting and much studied property of binary ionic mixtures, the linear mixing rule 
[M, IT2I IT3L IT^j. Hansen et al. have pointed out that the excess internal energy of the BIM, u = U ex /NkBT, is very 
acurately represented as a linear combination of the excess energies of two one-component plasmas with the same 
electronic coupling constant V as the BIM (and ionic charges Z\e and Z2&), i.e. the mixing of the two components 
at the same temperature and electronic density is largely ideal: 

u(T',x 1 )~x 1 uocp(r'zf)+x 2 uocp{r'zl) (13) 

We find that this rule is satisfied at very large asymmetries as well, with the largest relative deviations occurring 
at r' <C 1, i.e. small Au concentrations, in agreement with 0. Given this nearly ideal mixing behavior we assume 
that other system properties, e.g. viscosity, are bracketed by the component values as well. For an interpolation 
relation between the viscosities of the two OCP's at a given composition we now borrow some concepts from the 
linear transport theory of composite media. A common situation encountered in such systems is that of "impurities" 
with a generalized conductivity a\ and total volume fraction (p\ randomly dispersed in a matrix a 2 - Under these 
circumstances the effective medium theory |27| employs a mean-field like, self-consistent approximation to predict the 
medium conductivity on scales much larger than those of the inhomogeneities. For the case of viscosity this yields 
[2l| for the effective medium viscosity r\ m : 

£*-*Z*2l = o (14) 

We now note that the BIM is obtained by combining one-component plasmas with volume fractions <pi = ZiXi/(Z) 
and use the above relation to predict the mixture viscosity, essentially assuming that the theory applies for atomically 
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sized "impurity" domains. Using Eq. II II for the individual OCP viscosities 771 and 772 we obtain the results shown in 
Fig. El 

The degree of agreement with the simulation results is fairly remarkable. This may lead one to believe that the 
system is perhaps thermodynamically unstable and separating into two OCP phases 8J. However, we find no evidence 
for such a scenario and conclude that the decoupling signaled by the linear mixing rule along with the tremendous 
asymmetry between the ions lead to behavior mimicking that of a macroscopically mixed system. 

We now turn to the case of ionic interdiffusion in the BIM. Mutual diffusion in plasma mixtures plays an important 
role in the prediction of stellar structure [29j , as well as the stability of ICF targets [5j . For the case of a binary mixture 
the mutual diffusion coefficient D12 can be calculated in terms of the fluctuations of the microscopic interdiffusion 
current l30l: 



D12 = X1X2 



d 2 {f3G/N) 



dx\ 



D°u (15) 



P,T 



D° 12 = ^ / (ic(*)jc(0))dt (16) 
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where G is the Gibbs free energy. The thermodynamic prefactor that multiplies the Green-Kubo component D\ 2 
reduces to unity for dilute gas mixtures [3(|, but in low density, weakly coupled plasmas goes to (Z 2 }/(Z) 2 , which 
has been interpreted as an effect of the ambipolar electric field of the electrons [2^| . It is worth noting that the above 
relation is a good estimate for X1X2 [d 2 (f3G /N) / dx 2 ] p T at weak as well as strong couplings. Since the linear mixing 
rule holds well for all couplings, the change in Helmholtz free energy upon mixing at constant electronic density (and 
temperature) : 

AF s fi 

— = f(r', Xl ) - x^ocpP'Z?) - x 2 focp(F'Zi) (18) 
is very well approximated by the ideal entropy of mixing (with negative sign): 

(19) 
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If we assume that the system pressure p = Peiectronic + Pionic is entirely determined by the electronic density, i.e. 
Peiectronic 3> Pionic, which is consistent with the initial assumption of a rigid electronic background, r s <l, then AG = 
AF We can therefore immediately calculate the thermodynamic prefactor as (Z 2 )/(Z) 2 . The difference between 
BIM and ideal gas mixtures appears here to be an entropic effect induced by the charge neutralizing background, 
as mixing occurs at constant electronic density (i.e. constant electronic pressure) as opposed to constant molecular 
density (i.e. constant ideal gas pressure). 

In the course of the molecular dynamics simulations with various D + — Au 39+ mixtures we also calculated the 
microscopic interdiffusion current j c , and therefore were able to determine the Green-Kubo integrand D® 2 - The 
results are shown in Fig. 0] relative to Dq = a 2 ui p , along with the discussed thermodynamic prefactor estimate. We 
find that D® 2 is almost concentration independent for xau > 0.1, but appears to decrease fairly steeply at lower 
concentrations. The prefactor (Z 2 )/(Z) 2 has a simple behavior, with a sharp maximum for small amounts of highly 
charged ions [3l| . 

For a BIM with small ionic asymmetry kinetic theory estimates of the D® 2 were found to be in good agreement 
with simulations |26| . A simpler model, widely employed for astrophysics problems, was proposed by Paquette et al. 
[3^| . Its main assumption is that the Boltzmann equation can be used to calculate the transport coefficients of plasma 
mixtures modeled as BIM by making use of the Chapman- Enskog solution method j24[. The authors further argue 
that screened potentials, ZiZje 2 exp(—r/\)/r, are better suited for such estimates than the pure Coulomb interaction. 
In order to extend the validity of this approach to strong couplings they propose as appropriate screening distance A 
the larger of A73 and a, where Ad is the Debye screening length: 

knT 



X ° ~ keV + E^jJ (20) 
Under these assumptions the Boltzmann equation mutual diffusion coefficient is given in the first Enskog approximation 



as: 



Zk B T(Ah + M 2 ) 

[F>12\1 - T7TT- (il) 
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where f^ 1 ^ are collision integrals [24| that have been tabulated with high accuracy in Ref. (3^]. We perform such 
calculations for the D + — Au 39+ mixture using the slightly better second Enskog approximation |32|. The outcome, 
see Fig. 0J reproduces rather well the MD simulation results for D\ 2 at xau > 0.1, but not the mutual diffusion 
coefficient Z?i2- This is an important point that merits further discussion. In fact, there is no reason to expect that 
Chapman-Enskog estimates based on the Boltzmann equation for the ions can reproduce the full Dyi for either pure 
Coulomb interactions (with some reasonable cut-off) or screened potentials. As shown in j2|| for low density plasmas or 
more generally here, the prefactor value (Z 2 ) / (Z) 2 only arises when the electronic background is explicitly taken into 
account either through its ambipolar field in a dilute plasma kinetic description 26] or simply at the thermodynamic 
level in the context of the linear mixing rule. No such effect is included when the standard Boltzmann equation, which 
is consistent with ideal gas thermodynamics, is used to model the dynamics of the ions. Therefore it is reasonable 
to expect that such approaches can only provide estimates of the Green-Kubo part, i.e. Di 2 , of the mutual diffusion 
constant, as already evidenced by our simulation results and perhaps not fully appreciated before. We note however 
that the prefactor value quoted above may be a good approximation only for r s — > 0. For real systems, e.g. those 
encountered in ICF or astrophysics problems, r. departs significantly from zero and the electronic pressure and 
polarization effects can lead to phase separation [MEL especially at high asymmetries. The thermodynamic quantity 
(d 2 (/3G / N) / dx 2 ) p^t is connected to the spinodal decomposition line of the plasma mixture |3^|, and therefore proper 
estimates require careful calculations. In particular, it was pointed out in }33| that it is not sufficient to consider a 
BIM with compressible but non-polarizable electronic background since a consistent treatment can only be achieved 
with an appropriate modeling of electronic screening. 

The apparent failure at low Au concentrations of the screened potentials method proposed by Paquette et al. may 
appear at first puzzling since as xau decreases so does the effective coupling constant T e f f , and therefore the accuracy 
of the theory should in principle only increase. We note however that this also requires > a, which for our system 
and chosen T means xa u no bigger than ss 10 -3 . In fact, although not easily seen in Fig. 0] the theoretical values drop 
sharply for such compositions to values close to the MD result at the lowest simulated Au concentration, xau = 0.03. 
For this composition the screening distance a, although larger than A^, appears to be too small. This is perhaps not 
unexpected given that due to the significant separation between the highly-charged ions they are primarily screened 
by the small ions and the electronic background, which requires distances significantly larger than a when the charge 
asymmetry is very large. 

In summary, we have investigated using molecular dynamics simulations the viscosity and mutual diffusion coeffi- 
cients of a plasma model known as the binary ionic mixture (BIM) when the asymmetry is very strong. We discuss in 
light of the MD results important shortcomings of widely used theoretical models at large asymmetries. For viscosity, 
an OCP-based "one-fluid" theory proves inadequate for quantitative predictions and we propose a new method based 
on the linear mixing rule. A commonly employed model of ionic interdiffusion based on the Boltzmann equation 
compares reasonable well with the simulation results, but we point out that it provides only part of the mutual 
diffusion coefficient. The missing thermodynamic piece may be particularly important for large asymmetries and low 
concentrations of highly-charged ions, situations often encountered in both ICF and astrophysics applications. Finally, 
calculations that take into account electronic polarization effects are currently under way and will be reported in a 
future publication. 
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the auspices of the U. S. Department of Energy by University of California Lawrence Livermore National Laboratory 
under Contract No. W-7405-Eng-48. 



[1] Y. Rosenfeld, E. Nardi, Z. Zinamon, Phys. Rev. Lett. 75, 2490 (1995). 
[2] J.G. Clerouin, M.H. Cherfi, G. Zerah, Europhys. Lett. 42, 37 (1998). 
[3] G. Salin, J.-M. Caillol, Phys. Rev. Lett. 88, 065002 (2002). 

[4] J. A. Paisner, J.D. Boyes, S.A. Kumpan, W.H. Lowdermilk, M.S. Sorem, Laser Focus World 30, 75 (1994). 
[5] H.F. Robey, Y. Zhou, A.C. Buckingham, P. Keiter, B.A. Remington, R.P. Drake, Phys. Plasmas 10, 614 (2003). 
[6] P. Amendt, J.D. Colvin, J.D. Ramshaw, H.F. Robey, O.L. Landen, Phys. Plasmas 10, 820 (2003). 
[7] J.L. Milovich, P. Amendt, M. Marinak, H. Robey, Phys. Plasmas 11, 1552 (2004). 
[8] J. P. Hansen, G.M. Torrie, P. Vieillefosse, Phys. Rev. A 16, 2153 (1977). 
[9] J. P. Hansen, IR. McDonald, P. Vieillefosse, Phys. Rev. A 20, 2590 (1979). 
[10] J. P. Hansen, Phys. Rev. A 8, 3096 (1973). 

[11] N.S. Schulz, C. Canizares, D. Huenemoerder, K. Tibbets, Astrophys. J. 595, 365 (2003). 
[12] S. Ogata, H. Iyetomi, S. Ichimaru, H.M. Van Horn, Phys. Rev. E 48, 1344 (1993). 
[13] H. DeWitt, W. Slattery, G. Chabrier, Physica B 228, 21 (1996). 

[14] H.E. DeWitt, W.L. Slattery, G.S. Stringfellow, in Strongly coupled plasmas, edited by S. Ichimaru, p. 635 (Elsevier Science, 



6 



New York, 1990). 
[15] E.E. Salpeter, Aust. J. Phys. 7, 353 (1954). 

[16] C. Pierleoni, G. Ciccotti, J. Phys.-Condens. Mat. 2, 1315 (1990). 

[17] B. Bernu, P. Vieillefosse, Phys. Rev. A 18, 2345 (1978). 

[18] S. Ichimaru, S. Tanaka, Phys. Rev. Lett. 56, 2815 (1986). 

[19] J. Wallenborn, M. Baus, Phys. Rev. A 18, 1737 (1978). 

[20] G. Salin, J.-M. Caillol, Phys. Plasmas 10, 1220 (2003). 

[21] See, e.g., M. Deserno, C. Holm, J. Chem. Phys. 109, 7678 (1998). 

[22] M.P. Allen, D.J. Tildesly, Computer Simulation of Liquids, (Clarendon Press, Oxford, 1987). 
[23] Z. Donko, B. Nyfri, Phys. Plasmas 7, 45 (2000). 

[24] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-uniform Gases (Cambridge University Press, Cambridge, 
England, 1970). 

[25] J. P. Hansen, F. Jolly, I.R. McDonald, Physica 132A, 472 (1985). 
[26] D.B. Boercker, E.L. Pollock, Phys. Rev. A 36, 1779 (1987). 
[27] R. Landauer, J. Appl. Phys. 23, 779 (1952). 

[28] See, e.g., F. Lequeux, A. Ajdari, Phys. Rev. E 63, 030502 (2001). 
[29] S. Vauclair, Space Sci. Rev. 84, 265 (1998). 

[30] J. -P. Hansen, I.R. McDonald, Theory of Simple Liquids, 2 nd edition, (Academic Press, London, 1986). 

[31] The qualitative behavior of the thermodynamic prefactor remains unchanged (but its values are even hig her) if deviations 

from the linear mixing rule are taken into account and estimated using the relation proposed in Ref. |12|| (the present 

charge assymmetry is however larger than the ones studied in [T^l. 
[32] C. Paquette, C. Pelletier, G. Fontaine, G. Michaud, Astrophys. J., Suppl. Ser. 61, 177 (1986). 
[33] D. Leger, C. Deutsch, Phys. Rev. A 37, 4916 (1988). 



7 



T 1 1 — I I I I | 1 1 1 1 — I I I I | 1 1 1 1 — I I I I | 1 1 1 1 — I I I I | 



\ 




I ■ * ■ * ■ * 

10" 1 10° 10 1 10 2 

r 

FIG. 1: OCP viscosity: present simulations (circles), fit of simulation results - Eq. 1111 (solid line), predictive relation of Ref. 
(dashed line), Bernu-Vieillefosse interpolation formula Fnl (dot-dashed line), Chapman-Enskog estimate using screened 
potentials (dotted line). 
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FIG. 3: Viscosity of the D+ - Au 39+ BIM at T = 0.05 as a function of composition: simulations (diamonds), one-component 
model of Ref. [3 (dot-dashed line), one-component model using Eq. II II for the OCP viscosity (dashed line), two-component 
model Eq. 1141 ('solid line). 
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FIG. 4: Mutual diffusion constant contribution Di 2 (see text) of the D + — Au 39+ BIM at F = 0.05 as a function of composition: 
simulations (circles), screened potential model of Ref. ,32| (solid line). Thermodynamic prefactor (Z 2 ) / (Z) 2 (dashed line). 



